##########################################################################################
## https://sunduanchen.github.io/Scissor/vignettes/Scissor_Tutorial.html
library(data.table)
library(optparse)
library(ggplot2)
library(dplyr)
library(patchwork)
library(Seurat)

##########################################################################################

option_list <- list(
    make_option(c("--single_cell_file"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    gene <- "GKN1"
    single_cell_file <- paste0("~/20220915_gastric_multiple/dna_combinePublic/public_ref/singleCell/njmu/epiall_nor_PCA_50_RE0.5.Rdata")
    out_path <- paste0("~/20220915_gastric_multiple/dna_combinePublic/images/singleCell_MUC6_vln/GKN/")

}

##########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

out_path <- opt$out_path
gene <- opt$gene
single_cell_file <- opt$single_cell_file

dir.create(out_path , recursive = T)

##########################################################################################

sc_dataset <- load(single_cell_file, verbose = F)
sc_dataset_all <- epiall_nor_PCA_50_RE0.5
##Idents函数定义要取的类别
##########################################################################################
Idents(sc_dataset_all) <- "patient"

## JZ762突变未纳入分析
## JZ740、JR004、JR009、JR039是MSI、
sc_dataset_all <- subset(sc_dataset_all , idents=c("JZ710","JZ725","JZ732","JZ740","JZ751","JR004","JR008","JR009","JR039"))
#sc_dataset_all <- subset(sc_dataset_all , idents=c("JZ710","JZ725","JZ732","JZ751","JR008"))

Idents(sc_dataset_all) <- "sample"   

for( class in c("IM" , "IGC" , "DGC") ){

        print(class)

        ## 提取特定病理类型的样本
        sc_dataset <- subset(sc_dataset_all , idents=c(class))
        sc_dataset <- FindNeighbors(sc_dataset, dims = 1:10)
        sc_dataset <- RunUMAP(object = sc_dataset, dims = 1:10)

        p1 <- FeaturePlot(object=sc_dataset,features=c("GKN1" , "GKN2"),cols=c('grey','red'),label=F,pt.size=0.3)
  
        p2 <- DimPlot(sc_dataset, reduction = 'umap', group.by = 'celltype' , pt.size = 1.2, order = c(2,1))

        #out_name <- paste0( out_path , "/Scissor_STAD_" , gene , "_mutation.",class,".sample.",mut_type,".pdf"  )
        #pdf(out_name)
        p3 <- DimPlot(sc_dataset, reduction = 'umap', group.by = 'orig.ident',  pt.size = 1.2, order = c(2,1))

        plot <- p1 + p2 + p3 + plot_layout(ncol = 2)
        out_name <- paste0( out_path , "/" , gene , "_" , class , ".pdf"  )
        ggsave(file=out_name,plot=plot,width=12,height=9)
}


sc_dataset_all <- subset(sc_dataset_all , idents=c("JZ732"))
#sc_dataset_all <- subset(sc_dataset_all , idents=c("JZ710","JZ725","JZ732","JZ751","JR008"))

Idents(sc_dataset_all) <- "sample"   

for( class in c("IM" , "DGC") ){

        print(class)

        ## 提取特定病理类型的样本
        sc_dataset <- subset(sc_dataset_all , idents=c(class))
        sc_dataset <- FindNeighbors(sc_dataset, dims = 1:10)
        sc_dataset <- RunUMAP(object = sc_dataset, dims = 1:10)

        p1 <- FeaturePlot(object=sc_dataset,features=c("GKN1" , "GKN2"),cols=c('grey','red'),label=F,pt.size=0.3)
  
        p2 <- DimPlot(sc_dataset, reduction = 'umap', group.by = 'celltype' , pt.size = 1.2, order = c(2,1))

        #out_name <- paste0( out_path , "/Scissor_STAD_" , gene , "_mutation.",class,".sample.",mut_type,".pdf"  )
        #pdf(out_name)
        p3 <- DimPlot(sc_dataset, reduction = 'umap', group.by = 'orig.ident',  pt.size = 1.2, order = c(2,1))

        plot <- p1 + p2 + p3 + plot_layout(ncol = 2)
        out_name <- paste0( out_path , "/" , gene , "_" , class , ".MUC6MutSample.pdf"  )
        ggsave(file=out_name,plot=plot,width=12,height=9)
}
